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BILAYER PLATES: MODEL REDUCTION, T-CONVERGENT FINITE 
ELEMENT APPROXIMATION AND DISCRETE GRADIENT FLOW 

SOREN BARTELS, ANDREA BONITO 1 *, AND RICARDO H. NOCHETTO* 


Abstract. The bending of bilayer plates is a mechanism which allows for large deformations 
via small externally induced lattice mismatches of the underlying materials. Its mathematical 
modeling, discussed herein, consists of a nonlinear fourth order problem with a pointwise isometry 
constraint. A discretization based on Kirchlioff quadrilaterals is devised and its T-convergence is 
proved. An iterative method that decreases the energy is proposed and its convergence to stationary 
configurations is investigated. Its performance, as well as reduced model capabilities, are explored 
via several insightful numerical experiments involving large (geometrically nonlinear) deformations. 


1. Introduction 

The derivation of dimensionally reduced mathematical models and their numerical treatment is a 
classical and challenging scientific branch within solid mechanics. Various models for describing the 
bending or membrane behavior of plates are available, either as linear models for the description 
of small displacements or nonlinear models when large deformations are considered; see [11]. The 
development of related numerical methods has mostly been concerned with the treatment of second 
order derivatives and avoiding various locking effects. The rigorous derivation of the geometrically 
nonlinear Kirchhoff model for the description of large bending deformations of plates from three- 
dimensional hyperelasticity in [14] has inspired various further results, e.g., discussing other energy 
regimes in [15], or the derivation of effective theories for prestressed multilayer materials in [21]. 
Bilayer plates consist of two films of different materials glued on top of each other. The materi¬ 
als react differently to thermal or electric stimuli, thereby changing their molecular lattices. This 
mismatch allows for the development of large deformations by simple heat or electric actuation. 
Classical applications of this effect include bimetal strips in thermostats, while modern applica¬ 
tions use thermally and electrically induced bending effect to produce nanorolls, microgrippers, 
and nano-tubes; see [5, 18, 19, 22, 23]. Preventing undesirable effects such as dog-ears forma¬ 
tion, as reported in [1, 24], motivates the mathematical prediction of bilayer bending patterns via 
numerical simulation. This requires having a model as simple as possible to be ameanable to numer¬ 
ical treatment and analysis, but sufficiently sophisticated to capture essential nonlinear geometric 
features associated with large bending deformations. 

A two-dimensional mathematical model for the bending behavior of bilayers has been rigorously 
derived from three-dimensional hyperelasticity in [21]. It consists of a nonconvex minimization prob¬ 
lem with nonlinear pointwise constraint. The energy functional involves second order derivatives 
of deformations associated with the second fundamental form of the mid-surface. The pointwise 
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constraint enforces deformations to be isometries, i.e., that length and angle relations remain un¬ 
changed by the deformation as in the case of the bending of a piece of paper. A related numerical 
method has been devised and analyzed for single layer plates in [3, 4]. 

It is our goal to develop a reliable numerical method for the practical computation of large bilayer 
bending deformations. Our contributions in this paper are: 

• To present a formal dimension reduction model allowing for various effects not covered in the 
corresponding rigorous analysis in [20]; 

• To propose a discretization of the mathematical model and prove its T-convergence as discretiza¬ 
tion parameters tend to zero; 

• To construct a gradient flow method to compute stationary configurations; 

• To carry out several numerical experiments to illustrate the performance of our numerical method 
and explore the nonlinear geometric effects captured by the mathematical model. 

In the remainder of this introduction we discuss the mathematical model and our main ideas to 
deal with the ensuing strong nonlinearities. 


Description of bilayer plates. We consider a geometrically nonlinear Kirchhoff plate model that 
allows for bending but not for stretching or shear. This selection is related to the choice of a 
particular energy scaling, namely that the elastic energy is proportional to the third power of the 
plate thickness t. Given a domain uj C M 2 that describes the middle surface of the plate, the model 
formally derived in Section 2 consists of minimizing the dimensionally reduced elastic energy 

(1-1) E[y} = \ [ \H + Z\ 2 - I f-y 

^ J LO J UJ 

within the set of isometries y : uj —> M 3 , i.e., mappings satisfying the identities 
(1-2) [Vy] T Vy = / 2 <*=*> d i yd j y = 5 ij , i,j = 1,2, 

in uj and with prescribed values y = y d and Vy = on the Dirichlet portion dpuj of the 
boundary duj. Hereafter, I^ denotes the identity matrix in for d = 2,3, and H stands for the 
second fundamental form of the surface 7 = y (uj) parametrized by y with unit normal v, 

Hij = v ■ didj-y, v = dyy x d 2 y. 

The symmetric matrix Z is given and can be viewed as a spontaneous curvature so that in the 
absence of body forces f the plate is already pre-stressed. Given the identity for isometries 

\H\ 2 = \D 2 y\ 2 , 


we rewrite the energy E[ y] in (1.1) as follows: 
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(1.3) E[ y] := \ [ \D 2 y\ 2 + ^ [ didjy ■ 

- Ju ij=1 Ju 


l diy d- 2 y 

M3iy| X l&yl 




• y- 


It is tempting to simplify this expression further because |9;y| = 1, i = 1, 2, for isometries. However, 
we refrain from doing so for stability purposes anticipating that the isometry constraint will be later 
relaxed numerically. In particular, the normalization will enable us to prove various bounds for 
the variational derivative of the energy. In fact, notice that the energy E[ y] is finite for y £ 
H 2 (uj ) 3 n W^(uj) 3 such that |<9jy| > 1 , i = 1,2, as well as f £ L 2 {uj) 3 and Z £ L 2 (uj) 2x2 . The 
condition |<9jy| > 1, i = 1,2, will play a crucial role throughout this paper. We encode boundary 
conditions and the isometry constraint in the set of admissible deformations 

(1.4) A:= { y £ H 2 (uj ) 3 : y| 9dU = y D , S/y\ doU = $ D , [Vy] T Vy = I 2 a.e. in w} 
and define the tangent space of A at a point y £ A via 

(1.5) E[y] := {w £ H 2 (uj) 3 : w\ dr)UJ = 0,Vw|^ = 0, [Vw] T Vy + [Vy] J Vw = 0 a.e. in uj). 
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Note that it is always possible to extend the functional E to H 1 ^ cj) 3 as follows: 

E[y] := +oo, yeH 1 (u) 3 \A. 


Minimizing movements. To find stationary points of E in A. we propose a gradient flow in H 2 (u), 
i.e. according to the fF 2 -scalar product: if s G (0, oo) is a pseudo-time, we formally seek a family 
y G L 2 (0,oo;A) with d s y(s) G .F[y(s)] for s G (0,oo) and 


D 2 (d s y{s)) : D 2 w = -5E[y(s),w] 


Vw G .F[y(s)]. 


The expression SE[ y,w] stands for the variational derivative of E at y G A in the direction 
w G J-"[y(s)], which we make explicit below. We note that upon taking w = d s y(s) G -F[y(s)], we 
obtain formally 

^-E[y(s)] = - j \D 2 (d s y(s)) | 2 < 0; 

the energy thus decreases along trajectories. The formal gradient flow is highly nonlinear and re¬ 
quires an appropriate interpretation. We adopt the concept of minimizing movements and consider 
an implicit first order time-discretization of the f^ 2 -gradient flow via successive minimization of 

y ^WD'^y - y k )\\h^) + E[y] 

in the set of all y G A to determine y k+1 . Our motivation for the use of the H 2 metric to define the 
gradient flow is threefold. First, it simplifies the implementation since it leads to the same system 
matrix as the main part of the bilinear form associated with the bending energy (1.3). Second, 
it is sufficiently strong to provide control over discrete time derivatives and discretization errors, 
such as the isometry constraint, without imposing severe step size restrictions. Third, it may be 
regarded as a damping term modeling the deceleration of a bilayer within a viscous fluid, which in 
turn gives some physical interpretation. Since the nonlinear isometry constraint is treated exactly, 
this nonconvex minimization problem is of limited practical value. We thus propose instead a 
linearization of the isometry constraint which yields a practical scheme. 


Algorithm 1 (minimizing movement). Let r > 0 be the time-step size and set k = 0. Choose 
y° G A. Compute v fc+1 G -T[y fc ] which is minimal for the functional 

v ^ ^11 + E[y k + rv], 

set y k+1 = y k + rv fc+1 , increase k k + 1 and repeat. 

The linearized isometry condition included in the set J-”[y fc ] implies that for every admissible vector 
field v G J-”[y fc ] we have that the constraint residual of the corresponding update y k + rv satisfies 

[V(y fc + rv)] T [V(y fc + rv)] - J 2 = [Vy fc ] T [Vy A: ] - h + r 2 [Vv] T [Vv] > [Vy fc ] T [Vy fc ] - I 2 , 


where the inequality A > B for square symmetric matrices A, B means that A — B is semi-positive 
definite; in particular the diagonal entries of A and B satisfy an > bn for all i. Applying this 
formula inductively with [Vy°] T [Vy°] = I 2 , we see that [Vy fc+1 "] T Vy fc+1 > 1-2 whence \djy k+1 \ > 1, 
j = 1,2. Therefore, iffy^ -1-1 ] is well defined and the minimization problem in Algorithm 1 admits 
a minimizer. Moreover, the space -T[y fc+1 ] is well defined, even though y fc+1 may not belong to 
A. Hence, the iteration of Algorithm 1 can be repeated. Every iteration of the algorithm requires 
computing a solution v fc+1 G -T[y fc ] of the following Euler-Lagrange equation 

r f D 2 v k+1 : D 2 w + SE[y k + rv fc+1 , w] = 0 VwGJ r [y ?l ]. 
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In terms of the new iterate y k+1 = y fc +rv fc+1 this is equivalent to the nonlinear system of equations 


( 1 . 6 ) 


1 f D 2 (y k+1 — y k ) : D 2 w + [ 

J UJ JU3 


D 2 y k+1 : D z w 


f a fl ( diy k+l d 2 y k+l \ v 

—j- 


/ a,s, y t+1 . { ' 
J U) ^ L 


+ E / d ^y k+1 

i,j =i^ 

+ E [ didjy k+1 


i,j=l ' 


8\W 

\diy k+1 \ 

f Oiy fc+1 _ 

l|aiy fc+1 | X L|3 2 y fc+1 | 


diy fc+1 (diy fc+1 • giw) i d 2 y k+1 } 

|3iy| 3 J X |d 2 y fc+1 |J 13 

- d 2 w d 2 y k+1 (d 2 y k+1 ■ d 2 w) 


|3 2 y fc+1 | 3 


}Zij = [ 

J J UJ 


= / f-w, 


for all w G J-’fy^]. We show in this paper how to discretize this system in space and present an 
iterative algorithm for its approximation. We study these algorithms and employ them to compute 
several equilibrium configurations. 


Oidline of the paper. The paper is organized as follows. In §2 we discuss a formal derivation of 
(1.1) from three dimensional hyperelasticity, following a suggestion of S. Conti. In §3 we introduce 
Kirchhoff quadrilaterals , which is a nonconforming finite element specially taylored for this appli¬ 
cation. It is an extension of the Kirchhoff triangles [8, 3], and possesses the degrees of freedom for 
the function and its gradient at the vertices of the underlying partition Th', this facilitates imposing 
the isometry constraint at the vertices. It turns out that the space of deformations as well as that 
of discrete gradients are subspaces of i7 1 (w), which is extremely convenient to discretize (1.6). We 
next introduce space discretizations Eh of E and derive in §4 their T-convergence to E in H 1 {uS). 
As a consequence, we deduce convergence properties of discrete almost absolute minimizers. Then, 
we study a practical iterative algorithm for the solution of the space discretization of (1.6). We 
show convergence of such an iterative scheme, thereby proving existence and uniqueness of our 
fully discrete problem. We also prove several important properties of this gradient flow, such as a 
precise control of the deviation from the isometry constraint (1.2). We conclude in §6 with insight¬ 
ful numerical experiments. In fact, we compute several configurations (such as cylinders, dog ears, 
corkscrews) that are observed in experiments with micro- and nano-devices [1, 24]. We empha¬ 
size that some of these effects, e.g. corkscrew shapes, are obtained with anisotropic spontaneous 
curvatures Z and are therefore outside the framework developed and analyzed in [20]. 


2. Dimension Reduction: Bilayer Plate Model 


We consider a plate ut := u x (— i/2, i/2) C M 3 of thickness t > 0 and whose middle surface is 
given by w C K 2 as illustrated in Figure 1 (left). The plate is clamped on the left edge 8du and 
free on the rest of the boundary, and its length perpendicular to dr>LO is L. The upper and lower 
layers are composed of materials with different molecular lattices, for instance differing by a factor 
8t > 0; this could be achieved by thermal or electric actuation in practice. To understand the 
natural scaling between t and St, we assume that the middle surface of the upper layer contracts 
to length L(1 — St) whereas the middle surface of the lower layer expands to length L(1 + St), so 
that the plate ut bends upwards as in Figure 1. Due to the clamped boundary condition along 
one side we imagine that for small deformations the lower and upper middle surfaces are given by 
R± = 0~ 1 L( 1 ± St), as depicted in Figure 1 (right). Since we aim at capturing bending effects in 
the limit t —> 0, we impose the condition 


lim —— = lim —— = 
t —K) R-\- t—>0 li 


0 

L 


4 


K, 



k > 0 being the curvature. Since t = R + — R- = n 1 ((1 + St) — (1 — St)), we deduce 

, , , t 2 

(2.1) lim — = -. 

t —>0 St K/ 

It is thus natural to impose a linear scaling between t and St which involves the curvature that is 
expected in the limit of vanishing thickness for a pure bending problem. 



R + = R- + t 


Figure 1. Two layers of thickness i/2 are stacked on each other and form the bi¬ 
layer plate. The undeformed middle surface is denoted by u and deforms into the 
surface 7 . 

We consider the energy density W : K 3x3 x wt -> 1 

(2.2) W(F,x) = ^\F T F-(I 3 ±StN(x.')) T (I 3 ±StN(x.'))\ 2 ±x 3 >0, 

where x = (x', x' 3 ), I 3 denotes the identity matrix in M 3x3 , and | A| stands for the norm associated to 
the Frobenius scalar product A : B := Ylij =1 = H (H T F). Hereafter St is a parameter only 

depending on the thickness t, describing the lattice mismatch of the two layers and satisfying St ~ t, 
whereas N : lo — > R 3x3 is a symmetric matrix which encodes inhomogeneities (dependence on x') 
and anisotropy (rectangular molecular lattice rather than cubic and preferred directions) of the 
underlying materials. Together St and IV(x') describe the pre-stressed bilayer {x e 07 : ±x 3 > 0}. 
When St = 0 the two materials composing the bilayers reduce to one, the reference configuration is 
stationary in the absence of a force ff and thus stress-free, and the energy density becomes 

(2.3) W(F) = ^ | F t F — / 3 1 2 . 

This function is asymptotically equivalent to the simplest energy density that obeys the principles 
of frame indifference and isotropy, namely W(F) = W(QFR) for all Q,R £ SO( 3) and, see [15], 

(2.4) W{F) « dist 2 [F, SO{ 3)), 

where dist is given by the Frobenius metric. To see the relation between (2.3) and (2.4) we argue as 
follows. Let F be close to 50(3), which is to say F = Fo+eFi with Fo G 50(3) and F\ perpendicular 
to the tangent space Tp 0 5O(3) to 50(3) at Fo and e< 1; we thus deduce dist 2 (F, 50(3)) = e 2 |Fi| 2 . 
The space Tp 0 SO( 3) can be written as 

7> 0 5O(3) = {Z: F%Z + Z t F 0 = O}; 

this follows by differentiation of the condition F^ Fo = I. Consequently, the normal space Np 0 SO( 3) 
to Tp 0 SO( 3) is 

N Fo SO( 3) = {Y : Fq Y — Y T F 0 = 0}, 
as can be easily seen because Fp 0 5O(3) © Np 0 SO(3 ) = M 3x3 and 

Z : Y = tr ( Z t Y) = tr ((Z T F 0 )(F 0 r F)) = -tr ((F T F 0 )(F 0 T Z)) = -tr (Y T Z) = -Z:Y, 
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whence Z : Y = 0 and the subspaces are orthogonal. Since Tj £ Np 0 SO(3) we infer that 

| F t F - I 3 \ 2 = |(F 0 + eF^iFo + eT\) - I 3 | 2 = \2eF[F 0 + e 2 if T\| 2 = 4e 2 |Fi| 2 + o(e 2 ), 

which shows the asserted relation between (2.3) and (2.4) for small e. 

We are interested in the bending regime of the bilayer, which corresponds to energies comparable 
to the third power of the plate thickness, cf. [15, 21]. To a deformation u : ut —> M 3 of the plate 
we thus associate the scaled hyperelastic energy 

(2.5) I t [u]=t ~ 3 f (w(Vu, x) — i t ■ u) dx. 

Jl*)t ' ' 

The function f* is a body force, whereas the energy density W is written in (2.2) and reads 

W(F,x) = ^|f t f-m| 2 , 

with symmetric matrices M = M(pt),N = N (x) £ R 3x3 given by 


M : = 


Mil 


M 


M\2 
12 M‘22 


T 


:= h ± 2 S t N + 5 2 N 2 


N: = 


Nu 
T 


m 


m 


n 


Here and in what follows alternating signs correspond to the upper and lower layers in which we 


have X 3 > 0 and X 3 < 0, respectively, i.e., ±x 3 > 0. Moreover, Nu £ 
n £ M is constant (to keep the formal discussion simple), and 

M u = I 2 ± 25 t Nu + $t(N\i + mm T ), 


2x2 • 


is symmetric, m £ 


M 12 = ±25tm + 5 2 (iViim + nm), 
M 22 = 1 ± 2 5 t n + S 2 (n 2 + |mj 2 ). 


To derive a dimensionally reduced model we assume that the actual deformation u of the plate, 
subject to boundary conditions and outer forces, has the form 


u(x', x 3 ) = y(x') + x 3 b(x / ) 

with a vector field b : ui —> M 3 that is perpendicular to the surface 7 parametrized by y, i.e., we 
have dty • b = 0 for * = 1,2. In other words, fibers orthogonal to the middle surface in the reference 
configuration remain normal to 7 and deform linearly. This special form of u is consistent with 
[14, 21] for energy densities with a vanishing bulk modulus. In general, a more general expansion 
including quadratic terms in x 3 has to be used. 

With this ansatz we find that Vu = [9,;u] 3 =1 £ R 3x3 can be written as 

Vu= [V'y,b]+x 3 [V'b,0], 


where V ; stands for the gradient with respect to xh and deduce that 


hln] 



^|(Vu) T Vu - M\ z - f t • u) dx 



(V'y) T (V / y) — Mu 

-M? 2 


M\2 

b| 2 - M 22 


+ x 3 


(V / b) T V'y + (V'y) T V'b 
b T (V'b) 


(V / b) T b 1 


+ 4 


(V'b) T V'b 

0 


— fV • u > dx. 


In order for this integral to be bounded as t —> 0 we need that the term |6| 2 — M 22 be at least of 
order t 2 . Since we have <5* ~ t, this is guaranteed if we enforce that 

|b| 2 — (1 ± 25tn) — 5 2 {n 2 + |m| 2 ) = |b| 2 — (1 ± dt.n) 2 — 5 2 |m| 2 = —<5 2 |m| 2 , 

i.e., we impose the constraint 

|b| = 1 ± 5tn ± x 3 > 0. 
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Since b(x') = /3(x / )i/(x / ), where i'(x') := |^y(x')x^y(x')| ^he un ^ normal to the surface 7 at 
y(x'), we obtain that /3(x') = 1 ± 5tn which is for simplicity assumed to be independent of x' G ui. 
This in turn implies that 

(V'b) T b = 0, V'b = (1 ± 5 t n)Vu. 

Recalling that the hrst and second fundamental forms of 7 are given by 

(2.6) G = (V'y) T V / y, H = -(W) T V'y, 

and introducing 

G t ■= i- 1 ((V / y) T (V / y) - M u ) = ^{G - I 2 T - ^(iV 2 , + mm T )), 

K t := (V'b) T V'b = (1 ± 5 t n) 2 (VV) T VV, 
we infer that 

rr i_ 1 f f 1 \tG t -2x 3 {l±5 t n)H + xlK t -M 12 1 2 

,[ l -"5 " f ‘ 

Retaining only the terms of order 5f or lower, because the higher order terms vanish in the limit 
f -7 0, we obtain 

I t [u] « ^ J |^-(f 2 |G t | 2 + 4x1(1 ±5 t n) 2 \H\ 2 + xj\K t \ 2 

—4tx 3 (l ± 5tn)Gt : H + 2 tx 3 Gt : AJ — 4x1(1 ± 5tn)H : K t + 8<5f |m| 2 ^ — • u| dx. 

To ensure that the first term in the integral remains bounded as t -7 0 we require that 

G = I 2 , 

that is the parametrization y of 7 is an isometry. Consequently, we obtain 

G t = ^f2t~ 1 5t.Nn - f% 2 P, P := IV 2 ! + mm T . 



Since the quantities Nu, P , H , A', m are independent of X 3 , we carry out the integration over x 3 G 
(—t/2,t/2) and deduce that 

rt/ 2 

/ t 2 |G i | 2 dx 3 = 4t5 2 |Ahi| 2 + t<5 4 |P| 2 , 

J-t/2 

r t /2 +3 

/ 4x 2 (l±M 2 |tf| 2 dx 3 = -(l + <5 2 n 2 )|P| 2 , 

J-t /2 4 

p/2 +5 

y i/2 X 4 |R'| 2 dx 3 = -|Af, 

r t/2 

/ —4fx 3 (l ± <5 t n)G t : H dx 3 = 2t 2 «5 t iV 11 : if + f 2 5 3 nP : AT, 

J-t/2 

rt/2 +3X2 

/ 2tx\Gt : A"dx 3 =--^P : K, 

J-t/2 6 


p/2 +4 x 

/ 4x1(1 ± (5in)P : ATdx 3 = —^nH : K, 

J-t/2 8 


8(5 2 |m| 2 dx3 = 8t5 2 |m| 2 . 
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It remains to examine the forcing term ft, which we assume to be of the form 

f t (x',x 3 ) = t 2 f(x / , X 3 ) 

to give a nontrivial limit. In fact, if we let 

1 f */ 2 1 f t / 2 

f ( x 0 : = 7 / f (x 7 , x 3 ) dx’3, g(x') := - / x 3 f(x',x 3 )dx 3 

1 J-t/2 1 J-t/2 

then the contribution to the energy due to the body force becomes 

J b • u dx = J (f(x') • y(x') + g(x') • b(x')) dx'. 

Inserting these expressions back into // [u] . setting 

A := lim — G M 
t—>0 t, 

and keeping only terms of order one in t, we readily obtain 

lim it [u] - Y 2 J (\H \ 2 + 6 AIVn : Hj dx' + A 2 J (jiVn| 2 + 2 |m| 2 j dx' - j f • y dx'. 

If we further denote 

(2.7) Z := 3AIVii 

and ignore the second integral, which is constant and so independent of the surface 7 , we see that 
the dimensionally reduced model is governed by the energy 

E[y] = Y 2 jj H + z \ 2 dx - J f • y dx', 

where the parametrization y : oj — > M 3 of the surface 7 is an isometry, namely it satisfies ( 1 . 2 ). 
We remark that the derivation of the dimensionally reduced model can be carried out rigorously 
in the sense of T-convergence for a large class of isotropic energy densities [21]. The only required 
assumptions are the cubic energy scaling (2.5) and the proportionality <5 f ~ t of (2.1). The quantity 
— Z acts as a spontaneous curvature for the bending energy ,E[y] and specifies properties of the 
bilayer material. If the material is homogeneous and isotropic, then Z = a/2 with a £ R; we 
refer to [20] for a discussion of the qualitative properties of minimizers. On the other hand, the 
material could possess inhomogeneities and anisotropies which are x'-dependent and are encoded 
in JVn(x'); we discuss some options together with numerical experiments in § 6 . We observe that 
the components n and m of N play no role in the reduced energy. 

We assume that the plate is subject to clamped boundary conditions on a portion of du 

y = Yd, Vy = on d D u, 

where yo ■ w —> M 3 , : uj —> R 3x2 are sufficiently smooth, and <1?£) = Vy_o is an isometry in 

cu, i.e. <I?£)(x') T <l>£)(x') = I 2 for x' 6 oj. The variational formulation of the reduced plate model 
consists of finding y G A, defined in (1.4), such that 

(2.8) E[y} = \ [ \H + Z \ 2 dx' — [ f-ydx'. 

is minimized, where H is the second fundamental form defined in (2.6) and Z the spontaneous 
curvature of (2.7). The new scaling ^ is immaterial and just set for convenience. Existence of 
solutions of the constrained minimization problem is a consequence of the direct method in the 
calculus of variations. 



3. Kirchhoff Elements on Quadrilaterals 

The fourth order nature of (2.8) and the pointwise constraint (1.2) on gradients in the bilayer 
bending problem reveal that a careful choice of finite element spaces for spatial discretization is 
mandatory. To avoid C ll -elements, which are natural in H 2 but difficult to implement, we employ a 
nonconforming method that introduces a discrete gradient operator and which allows us to impose 
the constraint (1.2) at the vertices of elements. The components of the discrete deformation y/, 
belong to an H 1 conforming finite element space W/, and its discrete gradients to another H 1 
conforming finite element space GThe degrees of freedom of our numerical method are the 
deformations and the deformation gradients at the nodes of the partition Th of u into rectangles 
which are the vertices of elements. 


Definition 3.1. For a conforming partition Th of uj into shape-regular, closed rectangles with 
vertices Afh and edges £h we define the midpoints of elements and edges, the diameters of elements, 
and the maximal meshsize by 


z E : = 


1 

4 


£ 

zeA/^n T 


Z E ■ = 


1 

2 


£ ■. 

z£j\T h nE 


hr '■= diam(T), 


h = max hr 
TeT h 


for all T &Th and all E G Eh- For every E G £h we let n E ,t E G M 2 be unit vectors such that n e is 
normal to E and t e is tangent to E. We denote by z l E ,z 2 E G A 4 n E the end-points of E so that 
E = conv{z^,z|}. 

The following definition modifies the well known Kirchhoff triangles [ 6 , 8 ] to quadrilaterals and is 
related to [7]. We let Q r (T) and P r (T) denote the set of polynomials on T of partial degree r on 
each variable and of total degree r, respectively. 

Definition 3.2. Let Th = {T} be a partition of u C M 2 into rectangles as in Defintion 3.1. 

(i) Discrete spaces: Define 

W h '■= {wh G C(ul) : lOhlr £ Q 3 (T) VT G Th, Twh continuous in Afh, 

Vw h (z E ) ■ n E = i(Vi%(4) + Vw h (z E )) ■ n B VE £ £ h , }, 

G h := {fi h £ C{u) 2 : ^ h \ T G Q 2 (T) 2 VT £ Th}- 

(ii) Interpolation operator. Let X 2 : H 2 (lo) 2 —> Gh be defined by 

T 2 fi)( z) = -0(z) for all z G Afh, 

TlT>(z E ) = tf{z E ) for all E G £ h , 

^V’(zt) = | f° r al1 T G 7 ' h - 

z£j\f h nT 


The operator X) 2 is also well-defined for discrete vector fields if G VW. 
(Hi) Discrete gradients: Define Vh ■ —> Gh by 

V h w :=Ift[Vu:]. 

The operator Th. is also well-defined for discrete functions w G W h- 


Since the set of vertices, midpoints of edges, and element midpoint is unisolvent for the polynomial 
space <Q >2 (T), the interpolation operator is well-defined. However, Z 2 differs from the canonical 
nodal interpolation operator because of the condition at the element midpoints. The latter is 
imposed for practical purposes but makes Z^ inexact over Q 2 (T). Nevertheless, Z| is exact over 
Qi(T) so that the Bramble-Hilbert lemma implies 

(3.1) \\if - lhif \\ L p ( T ) + h T \\Tif - VZ 2 ^|| iP(T) < c 2 /i^||Z> 2 -i/:|| L p (r) 
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for all if E Wp(T ) 2 and 2 < p < oo. The operator Z| is also well-defined on VW/, since for 
every u>/, E W/, we have that V't/y, is continuous at the nodes Mh and at the midpoints of edges. 
We will also need the canonical nodal interpolation operator Z| : H 3 (uj) —> W h, which is defined 
by evaluating function values and derivatives at vertices of elements and normal derivatives at 
midpoint of edges by averaging. Since Z? is exact for w E P 2 (T), the Bramble-Hilbert lemma yields 

(3.2) ||w - Zfu;|| LP(T) + h T \\Vw - VZ^|| iP(T) + hff\\D 2 w - D 2 l%w\\ LP{T) < c 2 h^\\D‘ i w\\ LV{T) 
for all w E W 3 (T) and 2 < p < oo. A less obvious but useful stability bound reads 

(3.3) \\D 3 Zf L w\\i,p(T) < c\\D 3 w\\lp(t) for all w E W 3 (T) and 2 < p < oo. 

To see this, we first write Z() (w — q) = (w — «) + ( X h w — w^j for all (/ E P 2 (T). Therefore, invoking 
an inverse estimate together with D 3 q = 0, we use (3.2) to obtain 

II-D 3 Z^Hlp(t) < ch^ 2 \\S7Xl{w - <?)|| L p(t) 

< ch^ 2 \\V{w - q)\\ L p(T) + ch^ 2 \\V(liw - w)\\ LP{T) < c\\D 3 w\\ LP{T) , 

provided that q is appropriately chosen, e.g., as a suitable Lagrange interpolant of w over P 2 (Z). 
Hereafter, c > 0 indicates a generic geometric constant that may change at each occurrence, depends 
on mesh shape regularity, but is independent of the functions and parameters involved. 

Remark 3.1 (nodal degrees of freedom). The degrees of freedom in W/, are only the function values 
at the vertices (wh( z) : z E Mh), and the gradients at the vertices (Vid/,(z) : z E Mh)- In fact, the 
remaining four degrees of freedom of the finite element Q 3 (T) are the normal components X7wh(zE) 
of the gradients at the midpoints z^; of edges E which are fixed as the averages of directional 
derivatives \7wh{z l E ) -ng at the endpoints z l E of the edges. The values Vrc/j(z£;) -t e can be written 
in terms of Wh{z E ) and Vw/ l ( z l E ) t e for i = 1,2. The matrix realizing the operator : W h —*■ &h 
elementwise is required for the implementation of Kirchhoff elements. 



V ^ 



FIGURE 2. Schematic description of the discrete gradient operator V/,. Filled dots 
represent values of functions, circles of gradients, arrows of normal components, and boxes 
of vector fields. The normal derivatives in the cubic space on the left are eliminated via 
linearity. 


Remark 3.2 (subspaces of iL 1 (w)). Enforcing degrees of freedom of at vertices z E Mh f° r both 
function values and gradients implies global continuity; thus W h C iL 1 (w). Likewise, the degrees of 
freedom of G h at vertices and midpoints of edges guarantee global continuity; hence G h C H l {u) 3 . 

We collect important properties of the discrete gradient operator in the following proposition. 

Proposition 3.1 (properties of V/,). Let 2 < p < 00 . There are constants Ci, i = 1,..., 4, indepen¬ 
dent of h such that the following properties of the discrete gradient are valid: 

(i) For all Wh G W h we have 

(3-4) c^WVwhUp^ < \\VhWh\\LP(u) < Cl\\VWh\\LP(ujy, 
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(ii) For all Wh £ W h and, T &Th we have 

(3-5) c 2 1 \\D~ w h\\LP(T) < ||VV/ l w/ l ||i, P (7’) < C2||-D 2 w/ l ||LP('r); 

(in) For all iv £ Wp(T) and T £ Th we have 

(3.6) || V«; - Vft/wllip^) + h T \\D 2 w - \7V h w\\ LP ( T) < c 3 /4||D 3 «;|| LP ( T ); 
(iv) For all Wh £ W h and T £ Th we have 

(3.7) || Vw h - V h w h \\ LP{ T) < C4h T \\WhWh\\ L P(T)- 


Proof, (i) Given Wh £ W h the function = Vh^h £ G/j is well-defined and the operator : 
W h —>• G h is linear, whence = 0 implies F/h^h = 0. Conversely, if V h w h = 0 then we have 
that Vwh{ z) = 0 for all z £ A4 and Vu^(ze) = 0 for all E £ <5/,. Since the tangential derivatives 
of Wh vanish at the endpoints and midpoints of E £ £ h , and Wh is cubic on E, we deduce that Wh is 
constant on E. The fact that functions in Wh are globally continuous implies that Wh is constant 
over the skeleton of Th- Let T £ Th and note that there are four remaining degrees of freedom 
in <Q> 3 (T). Since Vwfcfzg) ■ n e = 0 for all E £ £h n T, we see that Wh is constant in T, whence 
Twh = 0. The equivalence of the identities Vwh = 0 and = 0 implies the asserted norm 

equivalence because Q:$(T) is finite dimensional. 

(ii) We proceed as in (i). If D 2 Wh = 0, then V't/y, is constant and so is V h^h according to its 

definition; thus VV hU>h = 0. Conversely, if = 0, then V hVOh is constant in T and thus Vu>/, 

is the same constant at the vertices and midpoints of edges of T. This matches the 16 degrees of 
freedom of Q 3 (T), whence Vwh is constant in T and D 2 iVh = 0. 

(iii) Estimate (3.6) follows from the interpolation estimate (3.1) with if = Tw upon noting that 

V hW . 

(iv) The estimate (3.7) is a consequence of (3.6), an inverse inequality, and (3.5). 

The independence of all constants of the element-size h t follows from scaling arguments. □ 

Remark 3.3 (bases of W h and G/J. We anticipate that our discrete algorithms (Algorithms 2 
and 3 below) do not require the choice of a particular basis for W h- Instead, we apply vertex 
based quadratures requiring only the values of the approximate deformation and its gradient at the 
vertices. In contrast, a basis for G h is required but standard. In our implementation we use the 
(nodal) Lagrange basis. 

4. Discrete Energies and T-Convergence of the Discretization 

We employ the Kirchhoff elements on quadrilaterals W? C iL 1 (w ) 3 and the discrete gradient op¬ 
erator Vh '■ W \ —>• G|, whose components are denoted d } ), j = 1,2, to approximate the energy E 
given by (1.3). For practical purposes, we also impose a relaxed isometry constraint at the vertices 
of elements, but we introduce a parameter 6 > 0 to control its violation. We will show in Section 5 
that in the context of an iL 2 -gradient flow, 6 is proportional to the gradient flow pseudo-timestep 
and can therefore be made arbitrary small. We next give a discrete version of (1.4) and (1.5). 

Definition 4.1. For 5 > 0, y o,h £ W 3 ( and &D,h £ ^h\d D n ^ the discrete admissible set be 

■A. 5 h := {y/i £ W l : yh\a D uj = yD,h\d D u>, V hyh\d D u = $D,h\d D u, 

|Vy ft (z)] T Vy ft (z) > I 2 Vz £ A f h , ||[Vy,] T Vy fe - / 2 || l i ( w) < 8}. 

The (pseudo) tangent space of A s h at y h £ A s h is defined by 
Th[yh\ := {w h £ W l : w h \d D co = 0, Vw h \d D w = 0, 

[Vw fc (z)] T Vy Jl .(z) + [Vy/i(z)] T Vwft,(z) = 0 Vz £ A f h ). 
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Notice that would be the tangent space to A s h at y if [Vy/ l (z)] :r Vy/ ) (z) were constant for 

vector fields in A'l at every node z G Afh, which explains the terminology. Notice also the use of 
the discrete norms which for 1 < p < oo are defined by 


\\< w : =Ej E ktwr. 

TeT h zeA^nT 

and satisfy the equivalence relation ||u/,J| i p( w ) ~ II^IIl^w) f° r piecewise bilinear functions Vh G 
C(cJ). We also define the discrete inner product (•, -)/i for piecewise continuous functions 4>, if G 

n Te 77 c°(T) 

0It(zMt(z) 

T&T h zeAT h rT 


and note that 1 )h- 

The finite element discretization Eh of the energy functional E in (1.3) is given by 


(4.1) E h [y h ] := - 


2 

|VV /l y fc | 2 + V (diim y h ] 


i,j=1 


d\ l yh d!]y h 
X 


|^y ft | \d$y, 


h | J 


,z t 




+ x(X Z)h — (f, y h)hi 

h Z 


for y/j G *4^ and Eh[yh] = oo otherwise, where T\ is the canonical Lagrange interpolation operator 
into the continuous piecewise Qi elements, and both Z and f are piecewise continuous in Q. The 
latter enables the use of quadrature for the last three terms, whereas the first term can be integrated 
exactly because V/,,y/,, is piecewise <Q> 2 - The energy (4.1) is thus practical. 


Remark 4.1 (discrete isometry relation). The nodal isometry relation [Vy/ 1 ,(z)] T Vy/ l (z) > I 2 
for yh G Ah implies that \dby h {x)\ > 1 for j = 1,2 and all z G Afh- Hence, the normalization 

djyh(z)/\djyh(z)\ in the discrete energy functional Eh[yh\ is well-defined. We will see that it 
allows for suitable energy bounds and gives rise to a coercivity property. 


We start by showing that the family {Eh}h >0 is (equi-)coercive. 

Proposition 4.1 (coercivity). Let the Dirichlet boundary data satisfy y £> G Lf 3 (w) 3 and G 
H 2 (u) 3x2 , and let yn.h ’■= Zf\yr>, &D,h ■= Let the data satisfy Z G (T ) 2x2 and 

f G n rerfe C' 0 (T) 3 . Let {y/),}/i>o be a sequence of displacements in H 1 (w ) 3 such that for a constant 
C independent of h there holds 

Eh[yh\ < C. 

Then yh G A s h and there exists a constant C depending on ||.Z’||l°°(u,); ||f and 

ll < ^ > D||j 7 2 (a;) > but independent of h, such that 

(4-2) ||VV h y h ||La (w) < C. 


Proof. We first argue that y h G A'l since otherwise Eh[yh\ = + 00 . As a consequence we have 
\dbyh(z)\ > 1 for j = 1,2 and all z G A 4, whence there exists a constant c independent of h such 
that 

(4.3) E h [ y h \ > -\\E r V h y h \\ 2 L 2 ^ - c(\\Vll[V h y h ] || L 2 ( w) \\Z\\ L oo^) + Hy^Hia^) ||f ||x,°°(a 7 ))- 

Since y h = y D,h and V^y^ = &D,h on dou, we can apply the Poincare inequality twice and bound 
11 y /t 11 l 2 (u) in terms of || VV^y/JI,^), ||yD,fe||#i( w ), and W^d^Wh 1 ^)- In view of (3.2) and (3.1), the 
latter two quantities are bounded by a constant times ||yD||// 3 (tj) and H^dIIh 2 ^), respectively. We 
observe that for all T G Th 

(4.4) |VXj[Vhyh]||x, 2 (T) < chT\\VlllV h y h \\\L°°(T) < chr\\VV h y h \\ La0 (T) < c\\X7V h y h \\ l?(T) 1 
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where the last step is an inverse inequality for £ G h- This implies the asserted bound. □ 


Remark 4.2 (coercivity and gradient flows). In the (energy decreasing) gradient flow setting 
adopted in Section 5, the assumptions of Proposition 4.1 are automatically satisfied provided the 
initial state has finite energy; see Proposition 5.2. 

We now show T-convergence of Eh to E in i/ 1 (w) 3 and deduce the accumulation of almost global 
minimizers of Eh at global minimizers of the continuous problem. For this, we assume that the 
discrete boundary conditions are obtained by interpolation of the continuous ones with strong 
convergence in L 2 (di)Uj). We also assume for simplicity that Z and f are piecewise constant. 

Theorem 4.1 (r-convergence). Let the Dirichlet boundary data satisfy yr> £ H 3 (u) 3 and £ 
H 2 (uj) 3x2 , and let yo,h '■= X , &D,h '■= X^d>£>. If Z and f are piecewise constant over the 
partition Th, then the following two properties hold: 

(i) Attainment. For all y £ A, there exists a sequence {y h}h with yh £ A® h C A s h for all h > 0 such 
that yh —» y in H 1 (u>) 3 and 

lim sup Eh[yh] < E[ y], 

(M)-> 0 

(ii) Lower bound property. Assume that 6 —> 0 as h —> 0. For all y £ iL 1 (w) 3 and all sequences 
{y/i} C H l (u) 3 such that yh —> y in H 1 (cu) 3 , we have 

E[ y] < lim inf Eh [y/i]. 

/i—>0 

Proof. We prove properties (i) and (ii) separately. 

(i) Since y £ A C H 2 (u) 3 , for every e > 0 the density of smooth isometries among isometries in 
H 2 (u ) 3 , cf. [17], implies the existence of an isometry y e £ H 3 (uj) 3 such that 

(4.5) ||y-ye||tf 2 (w ) < e. 

This, in conjunction with the isometry property of both y and y e , yields 

\E[y]-E[y e ]\<Ce. 

Therefore, we assume from now on that y £ H 3 (u) 3 and do not write the subscript e for simplicity. 
For h > 0 let y h = Z^y £ W'f be the nodal interpolant of y, i.e., we have y^(z) = y(z) and 
Vy^(z) = Vy(z) for all z £ A4- The latter yields [Vy/ l ] T Vy/ 1 = E at the nodes in A4, whence 
y h £ A^. Convergence of y h to y in H 1 {oj) directly follows from the interpolation estimate (3.2) 

lly - YhWmiu) < C2^ 2 ||y||//3 (a;) . 

It thus remains to prove the convergence of the discrete energies Eh[yh] to E[ y]. To derive the 
convergence of the first term in (4.1), we write 

(4-6) Vy - V h y h = V(y - X 3 y) + (VX 3 y - X2[VX 3 y]) 

and use (3.6) in conjunction with (3.3) to get 

||V(VX 3 y-X2[VX 3 y])|| L2(T) < ch T \\D 3 l 3 y|| i2(T) < ch T \\D 3 y || L2(T) . 

Combining this with (3.2) yields 

(4-7) [[VVfeyft, - D 2 y\\ L 2 (0J ) < ch||D 3 y|| L2M . 
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For the second term in Eh[ y/j], we first note that \d h -yh(x)\ = 1 for j = 1, 2 and all z G A4, and by 
nodal interpolation estimates 

• [d[ l y h x ^y fc ], Z*^ - {diX\[d^y h \ ■ [d^y h x d$y h ], Zy) 

< c X] /z T||- D (^[aj l y /l ])|| i 2 (r) ||D[^y h x d£y h ] \\ L 2 {T) 

( 4 - 8 ) TeTh 

+ cE ^T||^^[^jy/i]l|L 2 (T ) ||-D 2 [^fy/i X C>2yft]llL 2 (T), 

T£T h 

because D 2 diZf l [djyh\ = 0 for every T € Th and Z is piecewise constant over 7/,.; recall that 
c denotes a generic constant independent of h. Therefore, employing inverse estimates for both 
terms on the right-hand side of the preceding estimate, and recalling (4.4), we deduce 

{d&hi&jy h ] ■ [ d iyh x d 2 yh], z i:j ) h - (d^kidjy h ] • j^y .h x d$y h ], z^)| 

< c ^2 ^T||VVhyh||| 2 ( T )||Vhyft||i,<x>( T ) < ch\\V'Vhyh\\ 2 L 2( UJ )\\'Vhyh\\L°°{u)- 

T&T h 

We further observe that |j ^hyh\\L°°(uj) is bounded uniformly because y^ = Tf t y with y G H 3 (uj) 3 C 
C l (jU) being an isometry, and (3.4) with p = oo. This, together with (4.7), implies that the 
quadrature term above is bounded by ch\\y\\ 2 H3 ,y It thus remains to examine 

{dilhidj yh] ■ [o»i y ,h x d 2 y h ] , Zy) - (didjy ■ [dry x d 2 y \, Zy). 

Invoking again (4.6), we infer that ||Vy - V h y h \\ L 2 ^) < ch 2 \\D 3 y\\ L 2 ^ along with 

|| [d^y h x d%y h ] - [<9iy x d 2 y] ||l 2 (w) 

< II [0?y h - <9iy] X dfyyhWtffa) + II [^iy x [^y h - d 2 y] || L 2 (aj) < ch 2 ||D 3 y ||l 2 ( w ) 
because HV^y/Jl^oo^), || Vy || 

Z/0° (cj) < C. In addition, we see that 
VZhWhYh] - V 2 y = VZ 3 [V fe y h - Vy] + V(Z 3 [Vy] - Vy), 

along with 

||V(X 3 [Vy] - Vy)|| L2(T) < c^||X^ 3 y|U 2( ^ } . 

Using an inverse estimate and stability of in L°°(T), we get 

II VX^Vfcyh - Vy]|| L 2 (T) < ch^ l \\ll\W h y h - Vy]|| L 2 (T) 

< c\\ll[V h y h - Vy]||x,oo ( T) < c\\V h y h - Vy|| L oo (T) . 

Moreover, we further write 

W^hYk - Vy|| L oc (T) = ||X^[VXfy] - Vy|| L oo (T) 

< ~ VXjjy|| L ao (r) + ||V(X|y - y)|| L oc (T) , 

and obtain, according to (3.1) and (3.3) with p = oo and an inverse estimate, 

ll^lVX^y] - VXfy|| L oo (T) < ch T ||D 3 Xj]y|| £ 2 (T) < c/i T ||X> 3 y|| L 2 (T) . 

Since ||V(X, 3 y-y)|| L oo (T) < ch T \\D 3 y || L 2 (T) , we deduce ||VX^V^y/j]-V 2 y|| L 2 (T) < ch T \\D 3 y\\ L 2 {T) . 
This, together with the preceding bound for || [d±yh x d 2 yh] — \d\y x d 2 y\ \\l 2 (gj)i implies 

(diZlldfah] ■ [d^y h x d%y h ], Zy) - (didjy ■ [<9iy x d 2 y ], Zy) < c/?,||X> 3 y|| L 2 (aj) . 


(4.9) 
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Collecting the preceding estimates, we obtain [EyJy/j] — i£[y]| < ch\\D s y\\ L 2 / u \ where y is an 
abbreviation for y € . Selecting h = h(e) to be sufficiently small so that h\\D 3 y e \\ L 2 ^ < e yields 

\Eh[yh] ~ E[y}\ < ce. 
and the convergence of Eh[ y^] to E[ y] when /i -» 0 follows. 

(ii) We may assume that y^ £ ^ and Z/Jy/,] < C uniformly in h (perhaps for a subsequence not 
relabeled) for otherwise liminf/ l _ 5 .o Eh[yh] = +oo and there is nothing to prove. Hence Proposi¬ 
tion 4.1 implies that the sequence {Vhyh}h>o is uniformly bounded in i4(w) 3x2 . This guarantees 
the existence of 4? £ H 1 ( cu ) 3x2 such that after extraction of a subsequence (not relabeled) we have 
4>/, = VhYh — 1 $ in i4(w ) 3x2 and 4*/, -£ 4> in L 2 (u ;) 3x2 as h —> 0. The discrete approximation 
estimate (3.7) yields 


||Vy-$/»|| £ 2 (w) < ||Vy-Vyfe|| L 2 (aj) + ||Vy fc -< ||Vy-Vy/»|| £ a (w) +c/i||VVhyh|| £ 2 (w) , 

whence taking the limit when h 0 we deduce 4> = Vy and y £ H 2 (oj ) 3 because Vy h —»• Vy in 
L 2 (w ) 3x2 by assumption. Owing to the assumptions on the boundary data we have that y\d D u = yD 
and Vy|a DW = 4>d- To show that y is an isometry we utilize discrete interpolation estimates and 

y h € *4 


ll^h - ^2||L1( W ) < K *h - UiK ^]IUi( W ) + II Zhl*h $/J - ^2|| £ 1( W ) 

<cfe||V[$^ A ]|| £l(w) +co<5. 

The right-hand side converges to zero as (h, 5) —> 0 because of the uniform bound (4.2) of 4>/j in 
H 1 (cj) 3x2 . Hence, 4>^4>/j —>• I 2 pointwise almost everywhere in uj for an appropriate subsequence 
and, since 4>^4>^ —t 4> T 4> pointwise almost everywhere in uj, we deduce that y is an isometry a.e. 
in oj, i.e., y £ A. Since the H 1 -seminorm is weakly lower semicontinuous we get J' \D 2 y\ 2 = 
I ivif < liminf/,_ 5 .o f u |V4>/ 1 | 2 . It remains to prove that the following three terms tend to 0: 




r l ra.T 


1 ra.T. 


4 = \ diXh\djyh] 


diy h d%y h 

X 


, Zij 




13 ; > 


; > 


l\d?y h \ - |fl£y fc |J’“^ ~ ' [9 ' yh X ^ Zi 

II h = (diXl\d^y h \ ■ [d^y h x d%y h ], Z tj - {d^[d^y h \ ■ [d^y h x d%y h ], Z, 

III h = (dill[djy h \ ■ [d r {y h x d 2 y h \, - (d^y • [d x y x 3 2 y], Z, 

for all 1 < i,j < 2. We first note that 

141 < 11 d&h [®j"yh\ 11 £ 2 (w) 11 Zij 11 i°° (w)!! ld h yhl - \ d hy h \ - - - -qiL 2 ( W ) 

The first factor on the right-hand side is bounded as h —> 0 according to (4.4) and (4.2), and the 
second one by assumption. Since |<9,fy/iI > 1, we estimate the last factor as follows: 


B ' yi ' X - a}y h X 


sfy, x M -aty h x a 2 V* 


l^iYftl |^y/»| 


L K u ) 


< 


( d iyh ah ^ d Zyh . ah v / d \y h a \ 

fe “ d ' yh > x MM + 8iyh x ■ 82y ^ 


nh 

l^y^i 

3/). 


■I^y/i 




< l^yhl -1 £ 2 (w) + d lYh Li , u) \d 2 y h \ - 1 i4H . 


By the approximate isometry property and |clj l y/ l (z) > 1 for all z £ A 4, we obtain for y^ £ yl 1 



Moreover, since 'S/h.yh is uniformly bounded in w) 3 and w C l 2 , we have by Sobolev embeddings 
for all 1 < q < oo 

(4.11) ||V h yh|| L « (w) < c q \\V h y h \\ Lq{oj) < Cg(||Vhy ft || L 2 (w) + ||VV /l y /l || L 2 (aj) ) < c q , 

whence || \dby h \ — l|| L q^ < c q . Interpolating with discrete Holder inequalities between this discrete 
L q h -estimate and the discrete L^-estimate in (4.10), we deduce that 

IIi 8 M.i- 1 IIl»^ 0 . i = 1.2 

for p = 2,4 as <5 —> 0. This shows that \Ih\ — > 0. 

The second term Ilh accounts for the effect of quadrature and is the same as (4.8), whence 

\II h \ < c/i||VV ft yft|||2 (w) ||Vhyh|| L oo (w) . 

Since y h £ A S h is not an exact nodal isometry, we do not have direct control of 
IIV/j.y/j ||lo°(oj)- We invoke instead the two-dimensional discrete Sobolev inequality || V^y^H^o^) < 
c| log/i| 1 / 2 ||VV h y /t || L 2 M [10, p.123], to infer that 

\IIh\ < ch|log/r| 1 /2||vv,y,||| 2(w)/ ->0. 

The last term Illh is the same as (4.9) except that we do not have y £ H 3 (uj) 3 . We split Illh as 
follows: 

III h = ( {d l Xl [djy h ] - didjy ) • [<9iy x d 2 y \, Z tj ) 

+ ( d t T h[ d j l yh\ ■ { [d[ l y h X d%y h \ - [<9iy x d 2 y\ }, . 

We observe that dfLlfd^yh] —" didjy in L?{ cj) 3 , whence the first term tends to 0 as h —> 0. In 
fact, the uniform bound (4.4) on VX^[V/,y/,], in conjunction with (4.2), implies the asserted weak 
convergence, and the limit is found via 

(di(ll[d^y h \ - djy]),ip) = -{xl[d^y h \ - d^y h ],di^) - (d^y h - djy,diip) -A), 

which holds for every ip £ Hq(uj) 3 because 

WZhNhyh] ~ V ft y ft || L 2 (a;) < ch\\VV h y h \\ L 2 ^ < ch 

and V hyh —;> Vy in L 2 (uj). For the second term in Illh we resort again to the uniform L 2 -bound 
on VZ/jV/jy/J and write 

II [diyh x d%yh\ - [<9iy x d 2 y] IIl 2 ^) < ||V A y A - Vy|| L 4 (w) (||Vhyh|| L 4 (w) + ||Vy|| £ 4 (w) ). 

By compactness of the embedding H 1 (w) — y L 4 (cu), we have ||V^y^. — Vy||^4 ^ —» 0 as h —> 0. 
Finally, using (4.11) for q = 4 along with || Vy ||l°°(u;) X c because y is an isometry, we see that the 
preceding term tends to 0 and thus conclude the proof. □ 

Theorem 4.1 extends easily to piecewise constant approximations to L 2 -data Z and f and to piece- 
wise Lipschitz data over 77; we do not carry out the details. The following result is a consequence 
of standard abstract T-convergence theory [12, 9] combined with Theorem 4.1 and Proposition 4.1. 

Corollary 4.1 (convergence of absolute minimizers). Let 5 —> 0 as h ^ 0. Let C > 0 be a constant 
independent of h and {y/,}/, be a sequence of almost absolute discrete minimizers of Eh, namely 

(4.12) E h [y h \ < inf E h [w h ] + e h < C, 

w h&K 

16 



where e*. —> 0 as h —> 0. Then {y h]h is precompact in H 1 (w) 3 , and every cluster point y of yh is 
an absolute minimizer of E, namely 

(4.13) i£[y] = inf i?[w]. 

Moreover, there exists a subsequence of {yh}h (not relabeled) such that 

(4.14) lim ||y - = 0 and lim Eh[yh] = E[y\. 

h^f 0 h-> 0 

Proof. The uniform bound for the discrete energies and the coercivity property of Proposition 4.1 
imply that the sequence {Vy/,}* is precompact in L 2 (cj) 3x2 . Due to the norm equivalence (3.4) 
and a Poincare inequality we have that {y h}h is bounded in i7 1 (w) 3 . Moreover, because of the 
estimate (3.6) the differences V^y^ — Vyconverge strongly to zero in L 2 (w) 3x2 as h —> 0. Hence, 
there exists y € H 1 (w) 3 such that, up to the extraction of a subsequence, we have 

V/iy/i, Vyh —)■ Vy in L 2 (w). 

The lower bound assertion of Theorem 4.1 implies that y £ A and 

(4.15) E[y] < liminf E h [y h \. 

h-> 0 

It remains to show that y is a global minimizer of E. To prove this, let > 0 be arbitrary and 
z 6 A such that 

E[ z] < inf E[w] + p/2. 
we4 

The attainment property stated in Theorem 4.1 implies that there exist h > 0 and z h £ A s h so that 

Eh[zh\ <E[z\ + rj/2. 

On combining the previous two estimates and incorporating the fact that yh is a minimizer for Eh 
in up to the value £h, we have that 

Eh[yh] < E h [z h \ +£h< E[z\ + 1 ]/2 + £h < inf ^[w] + rj + e h . 

wGv4 

This together with (4.15) and the arbitrariness of p > 0 prove (4.14). □ 

Remark 4.3 (local minimizers). Statements about almost local minimizers of Ey, are not available 
in general. However, if E has an isolated local minimizer y, then there exist local minimizers {y/j.}/? 
of Eh converging to y provided h is sufficiently small [9, Theorem 5.1]. We defer the discussion of 
almost local discrete minimizers of Eh to Section 6. 

5. Fully Discrete Gradient Flow 

Corollary 4.1 guarantees that every accumulation point of almost absolute minimizers of {Eh}h 
is an absolute minimizer of E. We introduce and study in this section a practical gradient flow 
algorithm to minimize Ey,. on A') for h > 0 and where 5 is proportional to the gradient flow pseudo¬ 
time parameter. Our fully discrete gradient flow gives rise to an energy decreasing iterative scheme 
that converges to stationary points satisfying the isometry constraint up to a small error. However, 
like every gradient descent method, whether the algorithm reaches an almost absolute minimizer, a 
saddle point, or a local minimizer is not possible to discern. We discuss this further in Section 6. 

Algorithm 2 (discrete if 2 -gradient flow). Let r > 0 and set k = 0. Choose y° 6 Ay. 

(1) Compute y£ +1 6yj] + Eh [y|] which is minimal for the functionals 

y h ^||VV/ l (y/ l - y^)llz, 2 ^) + E h [y h ] 
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in the set of all y h G y£ + T h [y£]. 

(2) increase k —> k + 1 and continue with (1). 


Every step of the gradient flow requires solving a nonconvex minimization problem. Since the 
primary variables of interest are the discrete gradients 

■= V/,y^ +1 , $> h := VhYhi ■= V ft w h , 

with w h G T h [y£] , we let their columns be &h,ji ®h,j, ^h,j for j = 1 , 2 , and write the corresponding 
Euler-Lagrange equations as follows: 


(5.1) 


- <!>„], V* h ) + (V<& h , Vtt fc ) + J2 fah[*hj] ' [||^y x 


i Zij 


+ J2 fahl* h j] • x wh + ifE x = (f,w h ) h 

jj=l I ,21 | ,11 


for all Wft G Th [y*]. Hereafter, to have a simple and compact notation, we let P a be the operator 


Pa:= n( / 3 “rr? 

a V a 2 


for any given a G M 3 and observe that |P a | < 1 provided |a| > 1. Notice that we omit writing the 
time steps k and k + 1 in (5.1). Existence of a locally unique solution to (5.1) follows from a local 
contraction property of the fixed-point iteration defined in the next algorithm, in which we write 

t f M 

y h for y h ■ 


Algorithm 3 (fixed-point iteration). Let yh G Aff, define y 3 = y h, and set i = 0. 
(1) Compute <h £ +1 := V^y ^ 1 with y[ +1 G y h + T h [yh] such that 


(5.2) 


-(V[< +1 - $ h ], V* h ) + (V< + \ V* h ) = - E {dilhlVhj} 

i,j= 1 


l^,ll 


$ 


h, 2 


1 $ 


h,21 


> 


E (^[ $ : 


t 1 




^A .1 X 


K, 


h, 21 


$ 


h,l 


1 $ 


x 2 

r, o ’ 


h,l I 


,^ij) +(f,w h ) h 
/ h 


for all w h G P/Jy/,] with T h ■= Vm, 

(%) increase i -> £ + 1 and continue with (1). 


Note that (5.2) is a linear system for <1>^ +1 . Under a moderate condition on the step size r > 0 
the iterates are uniformly bounded and the map i —y is a contraction in a suitable P 1 -ball. 

In particular, the limiting discrete Euler-Lagrange equations (5.1) then admit a locally unique 
solution. This, and related properties, are discussed in the next two propositions. 


Proposition 5.1 (local contraction property). Let the mappings f, Z be elementwise continuous, 
let $ h = V h y h , and set C := max{l, ||V$h|| L 2 (a; )}. If 

B h ■■= {yh € Af? : ||VV h y h || L2(w) < 2 C}, 

then the nonlinear map >->■ < h |, +1 in (5.2) is well-defined from Bh into itself and is a contraction 
with constant | provided that r < Co, with Co > 0 depending explicitly on C, a Poincare constant 
cp > 1 of uo, ||f || an d H^Hl 00 ^)- Consequently, if Algorithm 3 is initialized with the k-th 
iterate of Algorithm 2, i.e. yh = y%, then the solution of Algorithm 3 converges to the unique 
solution of the Euler-Lagrange equation (5.1) within Bh. 
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Proof. We proceed in three steps and simplify the notation upon writing || • || = || • 

(i) Existence and isometry relation: Given y e h G Bh, and in particular y e h G Aff, we see that 
[& e h (z)] T & £ h (z) > 1-2 for all z G A4- We thus have |I > 1 for j = 1)2, so that the right-hand 
side of (5.2) is well-defined and the Lax-Milgram lemma gives the existence of a unique solution 
y £ h +1 Gy h + Fh [y/i] • Due to the definition of Fh [y/J we infer that 

H +1 (z) - $/ l (z)] T $ /l (z) + $h(z) T [$£ +1 (z) - $ fe (z)] = 0 
for all z G Mh- This implies y^ +1 G A™, namely 

K +1 (®)] T < +1 (z) = [S ft (z)] T * h (z) + [< +1 (z) - $,(z)] T [^ +1 (z) - <Mz)] > / 2 
for all z G A/),., because y^ G 71“. 

(ii) Uniform bound: We next show that the iterates satisfy ||V<h^|| < 2 C. For this, we choose 

*H = = V4y' +1 - yj 


in (5.2). This leads to 


I||vK +1 -$ h )f + |||V(^ +1 -<h ft )ll 2 + |||V ^ +1 || 2 - *I|V<M 2 

= (f,y ^ +1 - y h)h - Y ~ ®h,j\ ■ x h ' 2 


i,j =1 


M^lrl 


I K,2 1 


,3 




- £ (aji• [p 4 . - *M) x j#i + Sq x - *m)]. z ‘j) h - 

% yj = 1 I h j 21 I /i, 11 

Since |-P$«( Z )| < 1 for all z G A4, we deduce 

^l|vK +1 - $ ft )|| 2 < \\\^hW 2 + c(/,z )(|| y £ +1 - yj + ||vK +1 - $ fc )|| + l|v^||||< +1 - $ fc ||),. 


<!>■ 


with c(/, Z) > 0 only depending on f and Z. We now apply the Poincare inequality to ||y ^ +1 — y/i|| 
in conjunction with (3.4), namely ||y ^ ( +1 — ytill < cp||4 >^ +1 ~ where we let cp be the product 
of the Poincare constant of ui and the constant ci > 1 in (3.4) and take it to be cp > 1. Applying 
the Poincare inequality again, this time to T ;, +1 — we obtain 

I|| V (^ +1 - 4>,)|| 2 < ^IV^II 2 + c(/,Z) 4(2+ ||V$i||)||vK +1 - $ h )\\. 


To prove that ||V4 > ^ +1 || < 2(7 we assume by induction that ||V4>^|| < 2(7, which is valid for £ = 0. 
Using ||V$/i|| < C yields 

||VK +1 - $ h )|| 2 < tC 2 + 4r 2 c(/, Z) 2 c 4 P ( 1 + C) 2 . 


Since C > 1, choosing 


r < C\ := min j 


1 1 

2 ’ 16c(/, Z) 2 Cp 


< 


4c(/, Z) 2 c 4 


C 2 

pji + cy 


implies ||V(4 >^ +1 — 4>/i )|| 2 < 2tC 2 < C 2 , whence ||V4>^ +1 || < 2(7 as asserted. 
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(iii) Contraction property : It remains to show that the map i— y given in (5.2) is a contraction 

on Bh■ For this, we subtract the equations that define <P ^ +1 and in (5.2) and verify that 


-W +1 - & h ),v* h ) + (v($r 1 - n), v* h ) 


2 (f/ 


*j=i 


L Kii 








E (^L- 

i,j= 1 


$£ 2 $ftl 

■ M d>t d h,2 


-ft, 2 l 


-ft,ll 


, Zi 


+ e to* 


•hj j 


r *M ^ 

X 


*j=i 

2 


L $ 


!-ll 

ft,l I 


1 $ 


£-ii 

ft ,2 I 




*2 


+£ (adK:/] ■ be'i'M '.&•*«] • z 4' 


Bounding separately the sum of the first and third terms and that of the second and fourth terms, 
and using the admissible choice *ft = <P ^ +1 — $> e h , we deduce the estimate 

I||V(d>f 1 - ¥ h )\\ 2 < Ic 2 - 1 ||V « +1 - *£)||||V(*E - *i)\\, 

where we have used ||V*£|| < 2 C shown in step (ii), that |*^ ■(z)| > 1 for all f > 0 , j = 1 , 2 , and 
z G A/ft, and that for a, b G M 3 with |a|, |b| > 1 the following estimates hold 

Pa-Pb 

This implies the asserted contraction property with constant \ in Bh for t < Cq := min{Ci, C 2 }. 
Finally, if yft = y£ is the /c-th iterate of Algorithm 2, then any fixed point of (5.2) is a solution of 
(5.1) and conversely. This implies uniqueness of (5.1) within Bh- □ 


= 3|a — b 


a 

|a| 


b 

Ib| 


< la — b| 


Remark 5.1 (time step). If f = 0, then the preceding proof shows that r has to be sufficiently 
small so that r < (4c(Z) 2 c|>) , where c(Z) = \\Z\\ 

L°° (oj) and cp > 1 is the Poincare constant 
of cu with vanishing Dirichlet condition on dpu. Due to a repeated application of the Poincare 
inequality, the case f ^ 0 requires a stringent condition on r. 

The following proposition shows that the discrete H 2 gradient flow of Algorithm 2 is energy decreas¬ 
ing and becomes stationary, its iterates are uniformly bounded, and the violation of the isometry 
constraint is controlled by the pseudo-timestep size. 


Proposition 5.2 (properties of iterates). Let Z and f be piecewise continuous over the partition 
Th- Let {yft}fcl 0 C -A“ be iterates of Algorithm 2. We then have that for all k > 0 

k 

(5-3) E h [y k h +1 ] + ^ E II Wft(yf 1 - y£)|| 2 < E h [ y°], 

T e=o 


IIVVft,y^|| < C 
20 


and in particular 
(5.4) 



for a constant C > 0 depending on y°, / and Z, but independent of k. In addition, if 
[Vfty°(z)] T [Vfty°(z)] = I 2 for all z G A f h then y k h +l G A^ T , i.e. 

(5-5) ||[V h yE +1 ] T [V h y* +1 ] - / 2 || l i h < c 0 r Vfc > 0, 

where Co depends only on y°. 

Proof. The proof splits into three steps. 

(i) Energy decay: This is a direct consequence of the minimizing properties of the iterates, i.e., 

E h [y k h +1 ] + ^||VV,(y£ +1 - y^)||| 2M < E h [y k h }. 

(ii) Coercivity: For every y^ G Aff we have |<9j l y/i( z )| > 1 for 3 = 1,2 and all z G A4- This, in 
conjunction with (4.3) and (4.4), yields 

E h [y h \ > ^||VV h y h ||£ 2(h)) - c(\\Z\\ 2 LOO{u) + ||f|||oo M ). 

Since Eh[y k ] — Eh[ y°] from (i), this implies the asserted bound (5.4) of the iterates. 

(iii) Isometry violation: Abbreviating <F^ +1 := V/,,yJ[ +l and 4>| := V^y^, and noting that yj[ +1 — 
y k G E[ y|], we have for all z G A4, 

K +1 m - + [®{w] T K +1 w - = o. 

whence 

[*J +1 W] T l®S tl W] = + [(*E +1 - 4>k)(2)] T (($E +1 - 

A repeated application of this identity along with [V/jy°(,z)] T [V?jy°(z)] = I 2 for all 2 G My, yields 

IIK + 1 j T [ff; +1 ] - h Hm m < c^r nv k (yf+' - yi)ii 2 - 

e =0 

With the help of a Poincare inequality for V/ l (y^ +1 — yfj and (5.3), we deduce (5.5). □ 

We end this section by pointing out that the stationary state yff reached by the gradient flow 
might not be a discrete (almost) absolute minimizer. However, assuming that y?° is an (almost) 
absolute minimizer, then Corollary 4.1 guarantees that the accumulation points when h —>• 0 are 
absolute minimizers of the exact energy E. 

6 . Numerical Experiments: Performance and Model Exploration 

Algorithms 2 and 3 are implemented using the deal. II library [2], All systems of linear equations 
arising in the iteration of Algorithm 3 where solved directly using UMFPACK [13] leaving the discussion 
on developing efficient solvers open. The resulting deformations are visualized with paraview [16]. 
Note that only the displacement degrees of freedom at the vertices of Th are used for this purpose 
and the plates are thus displayed as continuous piecewise bi-linear elements. In addition, we say 
that a plate has reached numerically an equilibrium state parametrized by y^ +1 when 

(fU) febT 1 ] - gtbfll ^ 10 _e_ 

T 

where r is the gradient flow timestep. We set y“ := y^ +1 , where K is the smallest k that satisfies 
the above constraint. Each pseudo-time iteration k of the gradient ffow consists of a fixed-point 
iteration (Algorithm 3). This inner loop stops when the (£ + l)-th inner iterate satisfies 

)|VV ft (yf +1 -yf)||<<W P 
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for some <5 s t op > 0. It turns out that the number of subiterations experienced in practice is low (see 
for instance Figure 12). For ease of computation the discrete energies are approximated according 
to 

E h [y h ] ~\j \H h +Z\ 2 , 

where Hh is the approximation of the second fundamental form given by 

■= {{diYh)i x (dZyh) 2) • didjy h , 

without normalization of the vectors d^y^, i = 1, 2, and Z is given and symmetric. The absence of 
space dependence in Z corresponds to homogeneous materials. 

The purpose of the following numerical experiments is twofold. We first document the performance 
of Algorithms 2 and 3 by investigating their behavior for decreasing values of meshsize h and 
pseudo-timestep r and examining the violation of the isometry constraint. We also study various 
qualitative properties of the dimensionally reduced model, the existence of local discrete minimizers 
other than cylinders, which are for appropriate data and boundary conditions global minimizers 
according to [20], as well as the pseudo-evolution process (sometimes exhibiting self-intersections). 
We remark that our underlying nonlinear mathematical model allows for the description of large 
deformations which may be nonunique and cannot be expected to admit high regularity. Therefore, 
a meaningful error analysis seems out of reach. For small displacements the model leads to a 
linear bending problem for which best-approximation and interpolation results imply convergence 
rates. An experimental convergence analysis for the case of a simple exact solution reported below 
indicates a linear convergence rate. We document our quantitative findings in Tables 1, 2, and 
3 and use the symbol N/A to indicate special combinations of h and r for which we did not 
perform computations as these appeared to be irrelevant for the discussion and in some cases 
computationally expensive. 

Bilayer bending has technological applications in design and fabrication of micro-switches and 
micro-grippers as well as nano-tubes [5, 18, 19, 22, 23]. In these cases it is essential that the bilayer 
plate undergoes a complete folding to a cylinder without exhibiting dog-ears or a corkscrew shape, 
which may affect or impede the complete folding [24], Better understanding and control of this 
phenomenon is what motivated this work. We describe below several equilibrium configurations 
other than cylinders. 

6.1. Benchmark. We display in Figure 3 the pseudo-evolution of a bilayer plate oj = (—5,5) x 
(—2,2), clamped on the left-side dpuj = {x = —5} x [—2,2], i.e., 

y = 0, Vy = / 3x 2 on d D oj, 

with a spontaneous curvature Z = —12- The finite element partition consists of 5 uniform refine¬ 
ments of the rectangle oj. The parameters for Algorithms 2 and 3 are the pseudo-timestep r = 0.005 
and the sub-iteration stopping tolerance <5 s top = 10~ 4 . The discrete equilibrium state is a cylinder. 

6.2. Relaxation Process. We consider the clamped plate oj described in Section 6.1 for different 
spontaneous curvatures Z as well as different discretization parameters h. r. The subiterations 
stopping tolerance for Algorithm 3 is e> stop = 10 -3 . We examine the relaxation process towards 
equilibrium for two different spontaneous curvatures, namely Z = -1-2, —5/2- According to [20], 
the absolute energy minimizers are cylinders of height 4 (length of the clamped side) and of radius 
1 and 1/5 (reciprocal of the eigenvalues of Z) with an energy of 20 and 500, respectively. 

Table 1 documents the influence of meshsize h and pseudo-timestep r on the equilibrium shape 
and corresponding energy. The cases Z = —12 and Z = —5/2 are strikingly different. For Z = —12 
the equilibrium shape is a cylinder (absolute minimizer) provided r < CqH for a sufficiently small 
constant Co; see Figure 4 (left for r = 0.0025 and middle for r = 0.005). For Z = —5/2, and 
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FIGURE 3. Pseudo-evolution (clockwise) towards the equilibrium of a clamped rectangular 
plate with spontaneous curvature Z = —E- The bilayer plate is depicted (clockwise) for 
0.0, 0.1, 2.4, 60.0, 100.0, 130.0, 135.0, 207.3 times 10 3 iterations of Algorithm 2. The bilayer 
plate reaches a cylindrical shape asymptotically (limit of discrete gradient flow). This is an 
absolute minimizer. 


regardless of the size of r, the plate never reaches a cylinder but other equilibrium configurations 
(local minimizers) with much higher energy than 500; see Figure 4 (right). A plausible explanation 
is that the relatively large spontaneous curvature in the direction of the clamped side favors bending 
in such a direction, thereby creating a geometric obstruction to reaching a cylindrical shape. 

Table 1 also provides information about the threshold of r needed for convergence of the sub¬ 
iterations of Algorithm 3. This value, being sensitive to II-^IIl 00 ^), is more stringent for Z = —5/2- 
In fact, such iterations fail to converge for r = 0.02 and Z = —5/2, whereas for Z = —I 2 give rise to 
a local minimizer. This is consistent with Remark 4.3, which establishes the pessimistic thresholds 
To = 2.5-10 -3 for Z = —I2 and to = 10 -4 for Z = — 5/2 if we consider a Poincare constant cp = 10. 
In addition, we note that the case t = 0.00125 on the mesh ^6 and with Z = —I 2 yields a 
cylindrical equilibrium shape with energy of 17.2 (not reported in Table 1). This, in conjunction 
with the energies when t = 0.005 on the mesh resulting from four refinements and t = 0.0025 on 
the mesh resulting from five refinements, illustrates that Eh[ y“] — > E[ y] = 20 as /i —>■ 0, h = ct 
and c is sufficiently small to obtain cylindrical shapes. This is in accordance with Corollary 4.1. 


6.3. Asymptotics. In this section, we illustrate the predicted convergence rate of 0(t) for the 
violation of the isometry constraint in (5.5). We consider again the clamped plate uj = (—5,5) x 
(—2,2) described in Section 6.1 with a spontaneous curvature Z = — 12 . The space discretizations 
are subordinate to 4,5,6 and 7 uniform refinements of the initial partition of the plate consisting 
of 1 rectangle and are referred to as mesh $4, #5, #6 and #7, respectively. The sub-iterations 
stopping tolerance of Algorithm 3 is <5 s t op = 10~ 3 . The equilibrium isometry defect is defined to be 


ID h { y£°) := 


ll[V fe y^] T [Vhy^] 

M 
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^^\Mesh 

T 

#4 

4b 

II 

~h 

#6 

#7 

#4 

Z = 

#5 

—5/ 2 

#6 

#7 

0.02 

19.781n 

20.351n 

N/A 

N/A 

NoC 

NoC 

N/A 

N/A 

0.01 

19.335n 

20.157n 

20.576n 

19.590n 

575.372 

521.297 

536.036 

586.839 

0.005 

15.961y 

16.554y 

20.343n 

N/A 

567.886 

519.599 

534.365 

N/A 

0.0025 

15.765y 

16.395y 

17.304y 

18.062y 

581.405 

518.897 

N/A 

N/A 


TABLE 1. Equilibrium energies Eh[yf^] for spontaneous curvatures Z = —/ 2 ,— 5/2 and 
different meshsizes h and pseudo-timesteps r. The meshes correspond to 4, 5, 6, and 7 
uniform refinements of the plate to = (—5,5) x (—2,2). The symbol next to the energy 
values for Z = —J 2 indicates whether the equilibrium shape is a cylinder (y) or not (n) 
after the stopping test (6.1) is met. Typical equilibria are displayed in Figure 4 (left and 
middle, the latter corresponding to a local discrete minimizer). The numerical experiments 
indicate that the cylindrical shape (absolute minimizer) is reached for Z = —J 2 when r and 
h satisfy r < C^h for a sufficiently small constant Co- This experimental condition is more 
restrictive than the theoretically derived condition for mere convergence of our numerical 
scheme. The symbol NoC for Z = —5J 2 indicates that the sub-iterations of Algorithm 3 did 
not converge. The cylindrical shape is never reached for Z = — 5/ 2 ; see Figure 4 (right) for 
a typical equilibrium configuration. 



FIGURE 4. Equilibrium configurations of clamped rectangle plates u> = (—5,5) x (—2,2) 
described in Section 6.1 for different spontaneous curvatures and numerical parameters. Left: 
cylinder shape when Z = — J 2 , mesh refinement 6 and r = 0.0025; Middle: local minimum 
when Z = — J 2 , mesh refinement 6 and r = 0.005; Right: local minimizer when Z = —5J 2 , 
mesh refinement 7 and r = 0.01. If the spontaneous curvature is 1 or smaller, then the 
cylinder (absolute minimizer) is reached when the pseudo-timestep and the meshsize satisfy 
t < C 0 h for a sufficiently small constant Cq. For relatively large spontaneous curvatures, the 
curvature in the direction of the clamped side prevents the plate from bending completely in 
the orthogonal direction, thereby creating a geometric obstruction and leading to a (discrete) 
local minimizer for all numerical parameters tried. 


The results reported in Table 2 indicate that the predicted rate of convergence O(r) in Proposition 
5.2 is recovered numerically. This happens for both decreasing time steps r on a fixed mesh (columns 
1, 2, and 3) as well as simultaneous reduction of h and r while keeping h ~ r (diagonal). We stress 
that, according to the discussion of Section 6.2, not all equilibrium configurations are cylinders but 
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the experimentally observed linear rate applies to all of them. This is consistent with Proposition 
5.2 which is not specific to absolute minimizers. 


r 

#4 

#5 

#6 

#7 

0.02 

0.0559 

0.0488 

N/A 

N/A 

0.01 

0.0327 

0.0283 

0.0249 

0.0276 

0.005 

0.0180 

0.0155 

0.0139 

N/A 

0.0025 

0.0094 

0.0081 

0.0077 

0.0083 


TABLE 2. Isometry defect ID h { y£°) for the clamped plate u = (—5,5) x (—2,2) at equi¬ 
librium with spontaneous curvature Z = — 1^. A sequence of time steps r = 0.02 x 2 - *, 
i = 0,1,2,3, is considered for different space resolutions. A decay rate 0(t) is observed 
when the space discretization remains unchanged (columns 1, 2 and 3) as well as when the 
meslisize is reduced to satisfy h ~ t (diagonal). Notice that the isometry defect for a fixed 
time step is little affected by the space resolution (rows 2 and 4). 


For an experimental convergence test for the deformations we choose oj = (0, 27r) 2 , define 

1 0 

0 .5 


Z = — 


and consider clamped boundary conditions on the left side = {x = 0} x [0, 27t]. An exact 
solution is given by 

y (x,y) = (sin (x),y, 1 - cos(x)). 

Table 3 shows the scaled L 2 and H 1 errors for the stationary configurations computed with Algo¬ 
rithms 2 and 3; norms were computed with a one-point Gaussian quadrature rule on every element. 
The stopping criterion for the fixed-point iteration is <5 s t op = 10~ 4 . The underlying meshes corre¬ 
spond to £ = 3,4, 5,6 refinements of our coarse mesh so that hi = 27t 2 - ^. We choose a step size 
proportional to the mesh size, i.e., we set T£ = 2 _ ^/25 for £ = 3,4, 5, 6. The obtained errors indicate 
a nearly linear experimental convergence rate. 



#3 

#4 

#5 

#6 

IM/M 1/2 

IIVejIl/M '/ 2 

0.8250 

0.6723 

0.4273 

0.3622 

0.2310 

0.2002 

0.1220 

0.1077 


TABLE 3. Scaled approximation errors = y — y^° in an experimental convergence test 
with exact solution y given by a cylinder of radius 1. The approximations y“ are obtained 
with Algorithms 2 and 3 from a flat initial configuration and timestep sizes proportional to 
meshsizes. The numbers indicate a nearly linear experimental convergence rate. 


6.4. Effect of Aspect Ratio and Spontaneous Curvature. This section investigates numeri¬ 
cally the influence of (i) the spontaneous curvature (i.e. difference in material properties between 
the two plates) and (ii) the plates aspect ratio. The numerical parameters in Algorithm 3 are 
r = 0.005, (5 s top = 10~ 4 , and the partition corresponds to 5 uniform refinements of the initial plate. 
We consider plates u := (— L, L) x (—2, 2) with different lengths L > 0 and define the aspect ratio 
to be p := L/ 2. The plates are clamped on the side = {x = — L} x [—2, 2], It turns out that the 
tendency to bend in the clamped direction (accentuated for relatively large spontaneous curvatures 
in the clamped direction according to Section 6.3) is attenuated for small aspect ratios. We illustrate 
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this in Figure 5, which displays almost equilibrium configurations for p = 5/2,3/2,1,1/2 and 
spontaneous curvatures Z = — r /2 with r = 1,3,5. For large spontaneous curvatures, Algorithm 2 
did not always reach geometric equilibrium before the stopping test (6.1) was met. In addition, 
some pseudo-evolutions lead to severe folding and exhibit self-intersections, in which case they are 
no longer representative from the physics standpoint. 



16.55(20) 9.8(12) 6.5(8) 3.3(4) 


FIGURE 5. Equilibrium shapes of bilayer plates for several aspect-ratios p (from left to right 
p = 5/2,3/2,1,1/2) and spontaneous curvatures Z = —r /2 (from top to bottom r = 5, 3,1). 
Decreasing the aspect ratio restores the ability for the plate to fold into a cylindrical shape 
for larger spontaneous curvatures. For instance, this is the case for plates with parameters 
r = 3 and p = 3/2 or r = 5 and p = 1/2. Notice, however, that small regions around the free 
corners have not completely relaxed to equilibrium. This effect is due to the violation of the 
isometry constraint and reduced upon decreasing the discretization parameters as well as the 
stopping criteria. The numbers below each stationary configuration are the corresponding 
approximate energies. For comparison, the energies of corresponding plates with principal 
curvatures of - and 0 are given between parenthesis. 


6.5. Boundary Conditions and Shapes. We consider now different boundary conditions and 
plate shapes. We intend to examine the robustness of our numerical scheme in different situations 
and investigate plate shapes which are not studied in [20] and for which we do not know the absolute 
minimizers. 

Boundary conditions: We start with the plate co = (—3,3) x (—2, 2) clamped in a neighborhood 
of the bottom left corner, namely 8doj = {x = —3} x (—2,0) U (—3,0) X {y = —2}. We impose 
the spontaneous curvature Z = —I 2 , choose the numerical parameters r = 0.005, <5 s top = 10~ 4 in 
Algorithm 3, and use a partition of ui with 5 uniform refinements. Several intermediate shapes of 
the discrete gradient flow are depicted in Figure 6. The equilibrium configuration consists of a flat 
and a cylindrical part separated by a free boundary that connects points on the boundary at which 
boundary conditions change. 

Shapes: We now consider the I-shaped and the O-shaped plates depicted in Figure 7. The finite 
element meshes contain 7168 and 8192 quasi-uniform rectangles respectively. We set Z = — 5 / 2 , 
r = 0.005 and <f stop = 10 -3 . Relaxation toward numerical equilibrium shapes, which are not 
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FIGURE 6. Different snapshots of the deformed corner-clamped plate with spontaneous 
curvature Z = — I^. The equilibrium shape has energy 11.616 and is not a cylinder. It 
is worth comparing with the side-clamped plate discussed in Section 6.4, which leads to a 
cylinder of smaller approximate energy 9.81. 

cylinders, is depicted for both plates in Figure 8. We stress that for the I-shaped plate the curvature 
in the clamped direction dominates the other, thereby leading to a cigar shape that opens up at 
the bottom to accomodate the boundary condition. In contrast, the O-shaped plate is more rigid 
to bending and develops dog-ears at the free corners which prevent further bending. 
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FIGURE 7. Geometry of the I-shaped and O-shaped plates: the numbers indicate the length 
of the sides. In both cases, the clamped edge is the far left vertical edge. 


6 .6. Anisotropic Spontaneous Curvature. We now turn our attention to anisotropic sponta¬ 
neous curvatures, namely to matrices Z with different eigenvalues. In the first two examples the 
eigenvectors are aligned with the coordinate axis, but not in the third example. The spontaneous 
curvature is given by either 



with a = 1 or —5. The plate is u = (—2, 2) x (—3, 3) and the numerical parameters in Algorithm 3 
are r = 0.005, <5 s t op = 10 -3 . 

Dominant curvature: With a = 1 being the curvature in the clamped direction, we notice a rather 
minimal bending effect in such a direction. The plate pseudo-evolutions are displayed in Figure 9 
which shows an almost perfect rolling to a cylinder of energy 42.09. 

Curvatures with opposite signs: We take a = —5 to be the curvature in the clamped direction. This 
choice models the tendency of the plate to bend equally in opposite directions along the coordinate 
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FIGURE 8 . Different snapshots of the deformed I-shaped plate (left) and O-shaped plate 
(right). The corresponding stationary energies with spontaneous curvatures Z = —M 2 are 
404.57 and 314.152 respectively or about 12.448 and 14.137 relative to the plate areas. For 
comparison, the numerical stationary energy of a plate ui = (—5, 5) x (—2, 2) under the same 
boundary condition and spontaneous curvature is 519.6 or 12.99 once divided by the plate 
area (see Figure 5). It turns out that compared to the full plate, the stationary numerical 
energy per unit area is smaller for the I-shaped plate and greater for the O-shaped plate. 



FIGURE 9. Deformation of a plate with anisotropic curvature given by (6.2) with a = 1. 

The spontaneous curvature is 1 in the clamped direction, its effect being barely noticeable, 
whereas it is 5 in the orthogonal direction. The equilibrium shape is a cylinder (absolute 
minimizer) with an energy of 42.09. Compare with the case a = 5, presented in Section 6.4, 
for which the cylindrical shape is not achieved before the stopping test (6.1) is met. 

axes (principal directions). This is noticeable in the second and third snapshots in Figure 10, the 
latter also exhibiting self-crossing of the free corners. After three complete rotations, the plate 
relaxes to a cylindrical shape (absolute minimizer). Surprisingly, a cylindrical shape is reached 
before the stopping test (6.1) is met, unlike the case a = 5 (see first row - second column in 
Figure 5). 

Corkscrew shape: We consider now the second anisotropic spontaneous curvature Z in (6.2), which 
has eigenpairs 

Mi = 5 ei = [l,-l] T , M 2 = 1 e 2 = [l,l] T . 

This means that we still have principal curvatures 5 and 1 but with principal directions ei and e 2 
forming the angle 7r/4 with the coordinate axes. The deformation of this plate towards its equilib¬ 
rium shape is displayed in Figure 11. The plate exhibits a corkscrew shape before self-intersecting 
and continuing its deformation to a conic shape. In fact, a cylindrical shape is not reached before 
the stationarity test (6.1) is met. We emphasize that this is not in contradiction with [20], where 
scalar spontaneous curvatures are considered, and shed some light on equilibrium configurations 
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FIGURE 10. Deformation of a plate with anisotropic curvature given by (6.2) with a = —5. 

The spontaneous curvatures are —5 in the clamped direction and 5 in the perpendicular 
direction, which eventually dominates the former and leads to a cylindical shape after three 
full rotations. Snapshots are displayed counterclockwise, starting from the bottom right, for 
0.0, 0.3, 2.0, 4.0, 6.0, 10.0, 15.0, 18.0, 25.0, 172.0 times 10 3 iterations of Algorithm 2. The 
arrows indicate the clamped side. 

when the two principal spontaneous curvature directions are not parallel and orthogonal to the 
clamped side. Notice, however, that the equilibrium energy obtained is 61.31, which is larger than 
the cylindrical shape obtained when the principal direction aligned with the coordinate axes (see 
Figure 9). 




+■ 




FIGURE 11. Deformation of a plate with the second anisotropic curvature Z in (6.2). The 
principal curvatures are 5 and 1 but the principal directions form an angle 7t/4 with the 
coordinate axes. The snapshots are displayed clockwise starting at the top left, for 0.0, 0.1, 
0.4, 1.9, 3.0, 324.0 times 10 3 iterations of Algorithm 2. The plate adopts a corkscrew shape 
before self-intersecting. 


6.7. Energy Decay and Time Scales. One critical aspect missing in this study is the design 
of (pseudo)-time adaptive algorithms able to cope with the disparate time scales inherent to the 
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underlying energies. To illustrate this point, we plot in Figure 12 the energy decay of the clamped 
plate ui of Section 6.1 for spontaneous curvatures Z = —I 2 and Z = — 5 / 2 . Both energies exhibit a 
rapid decay at the very beginning of the deformation and very slow decay at the end. The numerical 
parameters of Algorithm 3 used for these simulations are r = 0.005, 5 s t 0 p = 10~ 4 and the finite 
element partition corresponds to 5 uniform refinements of the original plate. 



FIGURE 12. Energy decay versus pseudo-time for Z = —r /2 with r = 1 and r = 5. The 
cylindrical shape is reached when r = 1 after 207,000 pseudo-timesteps (total of 210.469 
solves accounting for the sub iterations). When r = 5, the equilibrium shape is reached faster 
after 86.600 pseudo-timesteps (total of 129.682 solves accounting for the sub iterations). 
However, the equilibrium reached is not a cylinder as already pointed out in Sections 6.2 
and 6.4; see for instance Figure 4. The energy decays fast at the very beginning of the 
relaxation process in both cases. In addition when r = 1, a second rapid decay arises with 
the unfolding in the clamped direction; see iteration 130,000 in Figure 3 (6th snapshot). 
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